{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding particles using NASA JPL Horizons system\n",
    "\n",
    "REBOUND can add particles to simulations by obtaining ephemerides from NASA's powerful HORIZONS database.  HORIZONS supports many different options, and we will certainly not try to cover everything here.  This is meant to serve as an introduction to the basics, beyond what's in [Churyumov-Gerasimenko.ipynb](../Churyumov-Gerasimenko).  If you catch any errors, or would either like to expand on this documentation or improve REBOUND's HORIZONS interface (`rebound/horizons.py`), please do fork the repository and send us a pull request.\n",
    "\n",
    "**Adding particles**\n",
    "\n",
    "When we add particles by passing a string, REBOUND queries the HORIZONS database and takes the first dataset HORIZONS offers.  For the Sun, moons, and small bodies, this will typically return the body itself.  For planets, it will return the barycenter of the system (for moonless planets like Venus it will say barycenter but there is no distinction).  If you want the planet specifically, you have to use, e.g., \"NAME=Pluto\" rather than \"Pluto\".  In all cases, REBOUND will print out the name of the HORIZONS entry it's using.\n",
    "\n",
    "You can also add bodies using their integer NAIF IDs:  [NAIF IDs](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/MATLAB/req/naif_ids.html).  Note that because of the number of small bodies (asteroids etc.) we have discovered, this convention only works for large objects.  For small bodies, instead use \"NAME=name\" (see the SMALL BODIES section in the [HORIZONS Documentation](https://ssd.jpl.nasa.gov/?horizons_doc))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'Sun'... \n",
      "Found: Sun (10) \n",
      "---------------------------------\n",
      "REBOUND version:     \t4.0.3\n",
      "REBOUND built on:    \tJan 12 2024 08:52:18\n",
      "Number of particles: \t1\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x1120e7440, m=0.9999999999950272 x=-0.007850021941778935 y=-0.003071795945987032 z=0.00020917793069980578 vx=0.00029283533072769824 vy=-0.0003990033211746936 vz=-2.86982605624576e-06>\n",
      "---------------------------------\n",
      "The following fields have non-default values:\n",
      "N:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1\u001b[0m\n",
      "python_unit_l:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 4097248214\u001b[0m\n",
      "python_unit_m:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2145773914\u001b[0m\n",
      "python_unit_t:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1864791206\u001b[0m\n",
      "rand_seed:\n",
      "\u001b[31m< 886170\u001b[0m\n",
      "---\n",
      "\u001b[32m> 676381\u001b[0m\n",
      "particles:\n",
      "\u001b[32m> (128 bytes, values not printed)\u001b[0m\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import rebound\n",
    "sim = rebound.Simulation()\n",
    "sim.add(\"Sun\")\n",
    "## Other examples:\n",
    "# sim.add(\"Venus\")\n",
    "# sim.add(\"399\")\n",
    "# sim.add(\"Europa\")\n",
    "# sim.add(\"NAME=Ida\")\n",
    "# sim.add(\"Pluto\")\n",
    "# sim.add(\"NAME=Pluto\")\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently, HORIZONS does not have any mass information for solar system bodies.  `rebound/horizons.py` has a hard-coded list provided by Jon Giorgini (10 May 2015) that includes the planets, their barycenters (total mass of planet plus moons), and the largest moons.  If REBOUND doesn't find the corresponding mass for an object from this list (like for the asteroid Ida below), it will print a warning message.  If you need the body's mass for your simulation, you can set it manually, e.g. (see [Units.ipynb](../Units) for an overview of using different units):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'NAME=Ida'... \n",
      "Found: 243 Ida (A884 SB) \n",
      "<rebound.particle.Particle object at 0x1120e7540, m=0.0 x=1.7283804595871177 y=-2.3772031668640636 z=-0.017147599198848423 vx=0.4522491120827006 vy=0.35453102610682247 vz=0.01092312049884117>\n",
      "<rebound.particle.Particle object at 0x1120e7540, m=2.1e-14 x=1.7283804595871177 y=-2.3772031668640636 z=-0.017147599198848423 vx=0.4522491120827006 vy=0.35453102610682247 vz=0.01092312049884117>\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dtamayo/Documents/workspace/rebound/rebound/horizons.py:167: RuntimeWarning: Warning: Mass cannot be retrieved from NASA HORIZONS. Set to 0.\n",
      "  warnings.warn(\"Warning: Mass cannot be retrieved from NASA HORIZONS. Set to 0.\", RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "sim.add(\"NAME=Ida\")\n",
    "print(sim.particles[-1])    # Ida before setting the mass\n",
    "sim.particles[-1].m = 2.1e-14 # Setting mass of Ida in Solar masses\n",
    "print(sim.particles[-1])    # Ida after setting the mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Time**\n",
    "\n",
    "By default, REBOUND queries HORIZONS for objects' current positions.  Specifically, it caches the current time the first time you call `rebound.add`, and gets the corresponding ephemeris.  All subsequent calls to `rebound.add` will then use that initial cached time to make sure you get a synchronized set of ephemerides.\n",
    "\n",
    "You can also explicitly pass REBOUND the time at which you would like the particles ephemerides:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'Venus'... \n",
      "Found: Venus Barycenter (299) (chosen from query 'Venus')\n",
      "Searching NASA Horizons for 'Venus'... \n",
      "Found: Venus Barycenter (299) (chosen from query 'Venus')\n",
      "---------------------------------\n",
      "REBOUND version:     \t4.0.3\n",
      "REBOUND built on:    \tJan 12 2024 08:52:18\n",
      "Number of particles: \t2\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x1120e7340, m=2.447838287784771e-06 x=-0.5245660051855766 y=-0.5089797719750951 z=0.023076250697312093 vx=0.8146329772804595 vy=-0.8461175012422552 vz=-0.05860561893434425>\n",
      "<rebound.particle.Particle object at 0x1120e75c0, m=2.447838287784771e-06 x=-0.6612400840226685 y=0.26970377066078294 z=0.04197305542931168 vx=-0.445669425641083 vy=-1.0953912971910784 vz=0.010725839790529726>\n",
      "---------------------------------\n",
      "The following fields have non-default values:\n",
      "N:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2\u001b[0m\n",
      "python_unit_l:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 4097248214\u001b[0m\n",
      "python_unit_m:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2145773914\u001b[0m\n",
      "python_unit_t:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1864791206\u001b[0m\n",
      "rand_seed:\n",
      "\u001b[31m< 829806\u001b[0m\n",
      "---\n",
      "\u001b[32m> 210256\u001b[0m\n",
      "particles:\n",
      "\u001b[32m> (256 bytes, values not printed)\u001b[0m\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "date = \"2005-06-30 15:24\" # You can also use Julian Days. For example:   date = \"JD2458327.500000\"\n",
    "sim.add(\"Venus\")\n",
    "sim.add(\"Venus\", date=date)\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the two Venus positions are different.  The first call cached the current time, but since the second call specified a date, it overrode the default.  Any time you pass a date (with time in UTC), it will overwrite the default cached time, so: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'Venus'... \n",
      "Found: Venus Barycenter (299) (chosen from query 'Venus')\n",
      "Searching NASA Horizons for 'Earth'... \n",
      "Found: Earth-Moon Barycenter (3) (chosen from query 'Earth')\n",
      "---------------------------------\n",
      "REBOUND version:     \t4.0.3\n",
      "REBOUND built on:    \tJan 12 2024 08:52:18\n",
      "Number of particles: \t2\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x1120e71c0, m=2.447838287784771e-06 x=-0.6612400840226685 y=0.26970377066078294 z=0.04197305542931168 vx=-0.445669425641083 vy=-1.0953912971910784 vz=0.010725839790529726>\n",
      "<rebound.particle.Particle object at 0x1120e72c0, m=3.0404326489511185e-06 x=-0.5575682509295975 y=0.813513379197337 z=0.00016705166681997484 vx=-0.8456273406927053 vy=-0.5626729140479645 vz=3.329436635861762e-05>\n",
      "---------------------------------\n",
      "The following fields have non-default values:\n",
      "N:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2\u001b[0m\n",
      "python_unit_l:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 4097248214\u001b[0m\n",
      "python_unit_m:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2145773914\u001b[0m\n",
      "python_unit_t:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1864791206\u001b[0m\n",
      "rand_seed:\n",
      "\u001b[31m< 433115\u001b[0m\n",
      "---\n",
      "\u001b[32m> 842290\u001b[0m\n",
      "particles:\n",
      "\u001b[32m> (256 bytes, values not printed)\u001b[0m\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "date = \"2005-06-30 15:24\"   # You can also use Julian Days. For example:   date = \"JD2458327.500000\"\n",
    "sim.add(\"Venus\", date=date)\n",
    "sim.add(\"Earth\")\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "would set up a simulation with Venus and Earth, all synchronized to 2005-06-30 15:24 UTC.  All dates should either be passed in the format Year-Month-Day Hour:Minute or in JDxxxxxxx.xxxx for a date in Julian Days.\n",
    "\n",
    "For reference HORIZONS interprets all times for ephemerides as [Coordinate (or Barycentric Dynamical) Time](https://en.wikipedia.org/wiki/Barycentric_Dynamical_Time)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Origin**\n",
    "\n",
    "REBOUND queries for particles' positions and velocities relative to the Solar System barycenter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'Sun'... \n",
      "Found: Sun (10) \n",
      "---------------------------------\n",
      "REBOUND version:     \t4.0.3\n",
      "REBOUND built on:    \tJan 12 2024 08:52:18\n",
      "Number of particles: \t1\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x1120e7e40, m=0.9999999999950272 x=-0.007850021941778935 y=-0.003071795945987032 z=0.00020917793069980578 vx=0.00029283533072769824 vy=-0.0003990033211746936 vz=-2.86982605624576e-06>\n",
      "---------------------------------\n",
      "The following fields have non-default values:\n",
      "N:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1\u001b[0m\n",
      "python_unit_l:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 4097248214\u001b[0m\n",
      "python_unit_m:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 2145773914\u001b[0m\n",
      "python_unit_t:\n",
      "\u001b[31m< 0\u001b[0m\n",
      "---\n",
      "\u001b[32m> 1864791206\u001b[0m\n",
      "rand_seed:\n",
      "\u001b[31m< 629898\u001b[0m\n",
      "---\n",
      "\u001b[32m> 447379\u001b[0m\n",
      "particles:\n",
      "\u001b[32m> (128 bytes, values not printed)\u001b[0m\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "sim.add(\"Sun\")\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reference plane is the ecliptic (Earth's orbital plane) of J2000 (Jan. 1st 2000 12:00 GMT), with the x-axis along the ascending node of the ecliptic and the Earth's mean equator (also at J2000).  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
